Computer simulation of topological defects around a colloidal particle or droplet 

dispersed in a nematic host 
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We use molecular dynamics to study the ordering of a nematic liquid crystal around a spherical 
particle or droplet. Homeotropic boundary conditions and strong anchoring create a hedgehog 
director configuration on the particle surface and in its vicinity; this topological defect is cancelled 
by nearby defect structures in the surrounding liquid crystal, so as to give a uniform director field 
at large distances. We observe three defect structures for different particle sizes: a quadrupolar one 
with a ring defect surrounding the particle in the equatorial plane; a dipolar one with a satellite 
defect at the north or south pole; and a transitional, non-equatorial, ring defect. These observations 
are broadly consistent with the predictions of the simplest elastic theory. By studying density 
and order-parameter maps, we are able to examine behaviour near the particle surface, and in the 
disclination core region, where the elastic theory is inapplicable. Despite the relatively small scale 
of the inhomogeneities in our systems, the simple theory gives reasonably accurate predictions of 
the variation of defect position with particle size. 

PACS numbers: PACS: 61.30.Cz,61.30.Jf,61.20.Ja,07.05.Tp 
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I. INTRODUCTION 

Suspensions of solid particles or immiscible liquid 
droplets in a host fluid (colloids) have attracted wide 
interest, due to a number of industrial applications (they 
appear in food, paints, ink, drugs) and fundamental re- 
search (they show features such as Casimir forces between 
the particles, formation of supermolecular structures and 
novel phases). Colloidal systems with a liquid crystal as 
a host fluid are of special interest [Q, |2|, [|. I n the ne- 
matic liquid crystals considered here, while the system 
remains translationally disordered and liquid-like, there 
is long-range ordering of molecular orientations about a 
preferred direction, termed the director. Elastic defor- 
mations of the director field around the colloidal parti- 
cles produce additional long-range forces between them. 
These interactions can be of dipolar or quadrupolar type 
depending on the symmetry of the director field around 
the particles [0 , and this in turn is extremely sensitive to 
the boundary conditions on the particle surface and the 
size of the particles j|, ^, |?| . 

To understand the interaction between the particles, 
it is important to know the liquid crystal ordering near 
one single particle. This ordering is also important for 
calculation of the elastic light scattering from nematic 
emulsions ||. The complication here is the presence of 
topological defects in such emulsions. Isolated drops pro- 
vide a spherical confining geometry for the liquid crystal. 
Sufficiently strong homeotropic anchoring of the director 
(that is, normal to the particle surface) then induces a 
radial hedgehog defect with topological charge +1. If the 
director field is uniform far from the particle, i.e. the total 
charge of the whole system is zero, topological consider- 
ations imply that an additional defect must be created 
in the medium to compensate the radial hedgehog. 

There are several types of defect which can arise in 
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FIG. 1: Sketches of the satellite defect and Saturn-ring de- 
fect. 



this case. Two are illustrated in Fig. [|. The first is a 
hyperbolic hedgehog with a topological charge —1, called 
a dipolar or satellite defect. The second is a quadrupolar 
or Saturn-ring defect, i.e. a —1/2 strength disclination 
ring that encircles the spherical particle. Theoretical and 
numerical work based on the elastic theory |)], [hJ suggest 
that the dipole configuration is stable for the micron- 
sized droplets which are usually realised experimentally; 
the Saturn-ring configuration should appear if the droplet 
size is reduced and, when present, it is always predicted 
to be most stable in the equatorial plane normal to the 
director. If the strength of the surface anchoring is weak 
enough, a third, surface-ring, director configuration is 
possible [|| |T^]. It has been also found that the twisted 
dipole configuration is possible for high enough ratios of 
twist to splay elastic constants JlO[ . 

However, the elastic theory has several disadvantages 
because of the way it treats the core structure of the 
defect embedded into the nematic host. First, the defect 
core region is considered as an isotropic inclusion with 
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some unknown free energy which is used as an adjustable 
parameter of the theory. Second, the liquid crystalline 
phase is treated as uniaxial. However, near the defect 
core (or indeed the particle surface) these assumptions 
are no lon ger valid. As has already been demonstrated 
fll| , |l2l , [l5, 14 , the defect core of the simplest disclination 
line of strength +1 has a rather complex structure: the 
order parameter varies near the core region; the uniaxial 
phase becomes biaxial; even the density of the mesophase 
can be affected by the presence of the wall or the defect 

For a better description it is necessary to use the full 
order parameter tensor to describe liquid crystalline or- 
dering [l?], |8). This tensor exhibits no anomalous 
behaviour at the singularity of the director field. More- 
over, one can resolve order parameter variation as well 
as biaxiality of the defect core using a Landau-de Gennes 
phenomenological expansion of the free energy. In this 
framework, only the liquid crystal density is assumed to 
be fixed. In its turn, density variation can be taken 
into account using Onsager theory fl5|j , although that 
approach will not be pursued here. 

Another approach to the study of defect core structure 
is to use computer simulation. Indeed, computer simula- 
tion is a well established method to treat bulk elastic co- 
efficients |y| ^o] , the surface anchoring strength jis], ^l) , 
and structures of disclination cores Jl4] , p2| . This means 
that computer simulation allows one to investigate the 
details of liquid crystalline structure which cannot be re- 
solved experimentally or using phenomenological theo- 
ries. 

In this paper we present the results of molecular dy- 
namics (MD) simulations of the topological defects which 
appear in the liquid crystalline mesophase around a 
spherical droplet or colloid particle. Our study comple- 
ments previous results obtained by Billeter and Pelcov- 



its |23| , who used systems of 2048 molecules modelled 
by the Gay-Berne potential. The present study inves- 
tigates the effect of changing the droplet size, and uses 
much larger numbers of molecules, up to one million. 
This turns out to be essential to avoid finite-size effects 
which are particularly important for the satellite defect, 
since it causes long-range director deformation. As well 
as observing both the Saturn-ring and satellite defects, 
with the latter being stable only for larger particle sizes 
as expected, we also observe a transition between these 
two structures involving an intermediate, non-equatorial, 
surface-ring defect. 

The paper is organized as follows. In Sec. || we present 
the computational details and molecular models we use to 
simulate the liquid crystal mesophase and the interaction 
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of the molecules with the droplet surface. Section 
contains the results of the simulations: density, director, 
order parameter, and biaxiality profiles for the Saturn 
ring, satellite, and surface ring defects. We also compare 
our results to the known predictions of the elastic theory. 
Concluding remarks are given in Sec. IV. 
II. MOLECULAR MODEL AND SIMULATION 
METHODS 



Molecular dynamics simulations were carried out using 
axially symmetric molecules interacting through the pair 
potential 
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(1) 



Here Qij = (rVj — cry + coVco; Ti i ^ s ^ ne cen ter-center 
separation, oq a size parameter, £0 & n energy parameter 
(both taken to be unity) and the orientation-dependent 
diameter er^ is defined by 



1 + X(u, -Uj) 



{iij ■ iij - ijj ■ Uj) 2 
1 - x(ui • Uj) 



-1/2 



(2) 



where \ = (k 2 — 1)/(k 2 + 1), n being the elongation. In 
this work we used k = 3 throughout. The orientation 
dependence is written in terms of the direction of the 
center-center vector fy = Vij/rij and the unit vectors 
iii, u.j which specify the molecular symmetry axes. This 
is a soft repulsive potential, describing (approximately) 
ellipsoidal molecules; it may be thought of as a variant of 
the standard Gay-Berne potential with exponents 

fi = 0, v = 0. The molecular mass m was taken to be 
unity, and the molecular moment of inertia fixed as / = 
2.5m<7g. 

The interaction of molecule i with the droplet was 
given by a shifted Lennard- Jones repulsion potential 



between the centers, having exactly the same form as 
eqn (|l|), but with replaced by Qi = (|r; — r c | — a c + 
<7o)/<To and a c = R + e>"o/2. Here, R is the colloid radius, 
and r c the position of the colloid center. This interaction 
potential results in homeotropic anchoring of the liquid 
crystal molecules, normal to the particle surface. 



The systems consisted of N — 8,000-1,000,000 par- 
ticles depending on the droplet radius (in the range 
R/ctq = 3 . . . 15) and defect type: details appear in Ta- 
ble |. A molecular dynamics program was run on a Cray 
T3E using a domain decomposition algorithm p3] . A re- 
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TABLE I: 


Run parameters 




R 


N 


Defect type 


Steps 


3 


8,000 


Saturn ring 


1,500,000 


5 


8,000 


Saturn ring 


1,500,000 


6 


8,000 


Saturn ring 


1,500,000 


7 


8,000 


Saturn ring 


1,500,000 


10 


64,000 


Saturn ring 


500,000 


10 


64,000 


off-center ring 


1,000,000 


15 


1,000,000 


Saturn ring 


500,000 


15 


1,000,000 


satellite 


500,000 



duced temperature k B T/e = 1 was used throughout (for 
this model, the phase behaviour is not sensitively depen- 
dent on temperature, as there are no attractive forces). 
The system size was chosen so that the number density 
of the liquid crystal far from the colloid was pa^ = 0.35, 
which lies well within the nematic region, with a corre- 
sponding bulk order parameter S ~ 0.81. For this sys- 
tem, in the reduced units defined by cto, £o and m, a 
timestep 8t — 0.004 was found suitable. 

Cubic periodic boundary conditions were employed in 
all the simulations. To facilitate the analysis, the posi- 
tion of the droplet was fixed in the centre of the simula- 
tion box, and a global constraint was applied to fix the 
orientation of the director along the z axis |27[] . 

To obtain the particular defect types of interest we 
have used several ways to prepare the initial configura- 
tion: 

1. The radius of the colloid droplet was steadily in- 
creased from a small value in systems which were 
already in the nematic state. This takes about 10 4 
steps to increase the droplet radius up to R/(Jq = 10 
and about 10 4 steps to equilibrate the system with 
the defect. This method generates an energetically 
favourable configuration of the kind which is stable 
for small colloid particles (i.e the ring defect). 

2. Disordered, isotropic configurations containing the 
colloid droplet were compressed to the ordered, ne- 
matic, state. This method, in principle, should give 
the lowest free-energy configurations in an unbiased 
way. However, for large systems, equilibration from 
the isotropic phase is extremely time-consuming: 
systems split into small domains which then realign 
(coarsen). 

3. For large systems, we prepared director configu- 
rations in the nematic state with the defect em- 
bedded in them. To introduce the defect into the 
mesophase we used approximate analytical expres- 
sions for the director field n(r) around the droplet, 



which were used to determine the initial molecu- 
lar orientations. This configuration was diluted, 
i.e. the molecule coordinates were scaled so that 
they did not overlap with each other and the col- 
loid droplet. We compressed the system under low 
temperature (T = 0.1) to the desired density. The 




FIG. 2: Coordinate systems used in this work. The system is 
assumed to have rotational symmetry about the z-axis. The 
director field is assumed to have axial (z) and radial (e p ) com- 
ponents only, and so may be described by the angle /3(r, 9) . 



low-temperature compression preserves the direc- 
tor structure in the system. Therefore, the original 
defect structure is still present in the system after 
the compression. After the compression, the sys- 
tem was equilibrated for 10 5 steps. This method 
allows one to prepare and observe metastable con- 
figurations, i.e. those which are local minima of the 
system free energy. 



III. SIMULATION RESULTS AND DISCUSSION 



Simulation results were analysed to give director, den- 
sity, order parameter, and biaxiality maps. Coordinate 
systems used to describe positional and orientational 
variables are shown in Fig. |2|. Positions, relative to the 
colloid particle at the origin, are expressed in spherical 
polar coordinates (r, 6*, 0); the system was assumed to 
be invariant to rotations about the z axis. Accordingly, 
all quantities were averaged over angle <f> (by rotating 
molecular positions and orientations together) before ac- 
cumulating simulation averages on a grid with Ar =0.05- 
0.1<7o, A9 = 10°. At each grid point, the order tensor Q 
was resolved into components based on cylindrical polar 
unit vectors, pointing (i) along the z-axis z, (ii) in the di- 
rection e p of increasing distance p = \/ x 1 + y 2 from the 
z-axis, and (iii) in the tangential direction = z x e p : 
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a,/3 = p, 0, z 



(3) 



where there are nujy molecules present in the bin {i,j}, 
S a f3 is the Kronecker delta, (...) denotes an ensemble 
average. Diagonalizing this tensor, for each bin, gives 
three eigenvalues and three corresponding eigenvectors. 
The largest eigenvalue defines the order parameter S(r, 6) 
for each bin. The biaxiality a is defined as 



a = g (S 2 - S 3 ) 



(4) 



where $2, S3 are the two remaining eigenvalues of the 
averaged order tensor Q [O]. 



A. Ring defect 

For all studied radii (R/ao — 3-15) the ring defect ap- 
pears immediately after equilibration of the system start- 
ing from the isotropic state or blowing up the droplet in 
the nematic state. We therefore concluded that this type 
of the defect is energetically more favourable for the cho- 
sen droplet sizes. For large configurations, we also used 
the method with the predefined director field. In the 
case of the ring defect, we used the trial function for the 
director angle (see Fig. ||) 



(3(r,9) = 



arccos \n 



■Xr,0)) 



proposed by Kuksenok et al || g|. This trial function 
captures the far-field behavior and a ring of (—1/2) sin- 
gularity at r — a r , 9 = ir/2: 



n, „n „ 1 sin 20 
0{r, 6) = 9 - - arctan 

2 



cos 29 



(5) 



Here f(r) is a unique function of r and independent of 
the angle 9. f(r) has also to comply with the boundary 
conditions: l//(oo) = (=> (3 = 0) far from the particle 
and f(R) = (=> (3 = 9) on the particle surface. The 
explicit form of the function f(r) can be found in refs 

A typical director map of the ring defect obtained from 
the simulations is shown in Fig. 0. From this map one 
can see that the ring defect does not have long-range di- 
rector distortions. Its core region is located very close to 
the droplet surface and the director distortion vanishes 
very quickly in the liquid crystal bulk. This is in agree- 
ment with the quadrupolar nature of the defect: far from 
the particle the director deviation angle has asymptotic 
behaviour (3 ~ (R/r) 3 sin 29. Therefore, this type of de- 
fect can be studied using a comparatively small number 
of particles and correspondingly small simulation boxes 
(see Table | for details). 




FIG. 3: Director map of the ring defect. Droplet radius 7? = 
15. The shading represents the value of the order parameter. 
The defect core is located close to the droplet surface, a r w 
1.164_R. The director variation rapidly vanishes in the cell 
bulk. 
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FIG. 4: Density profiles for the ring defect. We plot dimen- 
sionless number density pag versus reduced distance r/ao. 
Droplet radius R/ao = 15. The following directions are 
shown: 9 — 0, n (avoiding the defect) and 9 = n/2 (cross- 
ing the disclination ring). The inset shows the contour plot 
of the density map. 



Typical density profiles for the ring defect with R = 15 
are shown in Fig. ^ for 9 = 0, ir (not intersecting the 
defect) and 9 = it/ 2 (crossing the disclination ring). 
The profiles which avoid the disclination have oscillat- 
ing structure near the particle surface which is typical 
for a liquid-wall interface. The profiles which cross the 
disclination ring do not have oscillations. The difference 
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FIG. 5: Order parameter profiles for the ring defect along 
the directions 9 — 0,7r/2. Different curves correspond to the 
different droplet radii. The minimum of the order parameter 
defines the position of the disclination core. 



FIG. 7: Distance of the core region of the ring defect from 
the droplet center, as a function of reduced droplet radius 
R/ao. We also show the linear fit to the simulation results, 
a r = 1.164R - 0.33O-Q. 
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FIG. 6: Biaxiality profiles along the direction 9 — it/2 (across 
the defect). We plot the biaxiality parameter a versus re- 
duced distance r/ao. The peak of the biaxiality is centered 
on the defect core and coincides with the minimum of the 
order parameter (Fig ^|). 



may be due to the partial melting of the liquid crystal 
in the disclination core region. This melting damps the 
influence of the droplet surface on the interface region. 
Note that, even though the disclination core is located 
close to the particle surface, the region where the density 
has an oscillating structure is quite big. Therefore, the 
simulation box has to be of adequate size in order not to 
affect the defect structure. 

The order parameter profiles for 8 = 0, n/2 are shown 
in Fig |5[ The shape of these profiles in general reflects 
the typical structure of the core: the center of the core 
has lower order than the bulk and the core region extends 
over a few molecular lengths. Using order parameter pro- 
files one can also define the position of the disclination 
ring as the location of the minimum of the order param- 
eter. For large enough R a second minimum appears, 



indicating that the disclination ring defect has a com- 
plex structure. This minimum is probably due to the 
density modulation near the droplet surface. As we will 
see, the density modulation also influences the rest of the 
profiles. 

To emphasise the complex structure of the defect core 
we plot the biaxiality profiles (Fig ^j). From the order and 
biaxiality profiles one can see that the main biaxial ring 
is accompanied by one additional biaxial ring for large 
R, or even two extra rings for some R (e.g. R/ao = 5,6). 
Again, this is most likely due to the density modulation 
near the droplet surface. 

Using the minima of the order profiles (or, with the 
same results, the maxima of the biaxiality profiles) we 
extracted the distance from the core of the disclination 
ring to the particle surface. The dependence of this dis- 
tance on the particle radius R is shown in Fig [?]. It is 
interesting to compare with the phenomenological the- 
ory predicting linear dependence of the ring defect ra- 
dius on the droplet radius: a r rs 1.25i? from minimiza- 
tion of the elastic free energy using a trial function [p|, 
a r ?s 1.13i? using a simulated annealing method 
or a r w 1.26i? |l(J. Our MD simulation results give 
a r - R = -0.33cto + (0.164 ± 0.004)i? which is in good 
agreement with the phenomenological theory, especially 
if one bears in mind the complex structure of the defect 
core. 



B. Satellite defect 

The director field around the satellite defect is illus- 
trated in Fig [|. The director distortion extends much 
further than that of the ring defect. This is a direct con- 
sequence of the symmetry of the director distribution of 
the satellite defect: far from the particle, the director an- 
gle vanishes as j3 ~ (R/r) 2 sin 6, i.e. it is like a dipolar 
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FIG. 8: Director map of the satellite defect. Droplet radius 
R/ao = 15. The background contour plot represents the value 
of the order parameter. 




FIG. 9: Density profiles for the satellite defect. Droplet 
radius R/ao = 15. The following directions are shown: 6 — 
0, 7r/2 (both avoiding the defect) and 9 = ir (crossing the 
defect). The contour plot of the density map is shown in the 
inset. 



term in a multipole expansion. Therefore, to study the 
satellite defect, one needs very large systems. We used 
one million particles and droplet radius R = 15. 

The approach involving a trial function was used to 
generate the dipolar initial configuration, as in the case 
of the Saturn ring. In this case superposition of the radial 
and hyperbolic hedgehogs implies the following ansatz 



(3(r, 9) — 9 — arctan 



l//(r) + cos0" 



(6) 



Again, f(r) is supposed to be a unique function of r which 
does not depend on the angle 9. The boundary condi- 
tion now is l//(oo) = (=> (3 — 0) far from the parti- 
cle; the position of the hyperbolic hedgehog is defined by 
f(a s ) = 1. After the initial configuration was generated 



FIG. 10: Order parameter profiles for the satellite defect 
along the directions 8 = 0,7r/2,7r. Droplet radius R/ao — 15. 
The defect core is elongated with the long axis parallel to the 
far-field director. 
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FIG. 11: Biaxiality profiles for the satellite defect (9 — 
0, 7r/2,7r). Droplet radius R = 15. From these profiles one 
can define the center of the satellite defect as a maximum of 
the biaxiality. The small blip close to the droplet surface is 
due to the surface effects: it is seen on all profiles. 



the system was equilibrated for 50, 000 steps with the 
subsequent production runs extending to 500, 000 steps. 

The density profiles for the satellite defect are shown 
in Fig ^. The following directions are shown: 9 = tt 
(through the defect) and 9 — 0,7r/2 (avoiding defect). 
One can see that the density profile at 9 = n (across the 
defect core) has less prominent oscillating structure than 
the other two: this is again presumably due to the partial 
disordering of the mesophase in the core region. 

The order parameter profiles for 9 — 0,ir/2,ir are 
shown in Fig [l(]. One can see that the defect core re- 
gion is quite elongated with its long axis parallel to the 
far-field director. It is difficult to judge the position of 
the satellite defect from the order parameter profiles: the 
region occupied by the disclination core is rather large. 
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The biaxiality profiles for the same values of 9 are 
shown in Fig [ll]. From these profiles one can see that 
the center of the defect core is located at a distance 
a s « 1AR. The value predicted by the elastic theory 
is about \.22R for the simulated annealing calculations 
||; 1.17-R or IA&R §] for the free energy minimiza- 
tion using a trial function. Again, as in the case of the 
ring defect, both order parameter and biaxility profiles 
are affected by the density oscillations near the particle 
surface. 



C. Off-center ring 




Simulation results show that both the satellite and ring 
defects are at least metastable for R/cr = 15: once the 
particular defect is realized in the system, it is stable 
over the timescale that is accessible to our simulations. 
However, the satellite defect is not stable for the smaller 
droplets. Indeed, we observed a rapid transition (several 
thousands of MD steps) of the satellite defect to the ring 
defect for R/<tq < 10. Equilibrating the initial configu- 
ration with the satellite defect in the cell with droplet 
radius R/(Tq = 10, we observed that it evolves into an 
off- centered ring defect: a typical nematic director map 
is shown in Fig 12, The ring moved slowly, evolving to- 
wards an equatorial Saturn ring configuration. Doing 
long runs (up to a million timesteps) we conclude that 
this is an intermediate state between the satellite and the 
Saturn ring defect. The z-coordinate of the ring, relative 
to the droplet center, as a function of the number of the 
timesteps, is shown in Fig [H?]. The evolution dynamics is 
quite slow. In principle, one might observe the same type 
of transition in the opposite direction, for droplet sizes 
large enough to stabilize the satellite defect; this is still 
unreachable for the system sizes we explore here. Note 
that elastic theory also predicts the off-center ring con- 
figuration to be unstable ]To| , with the transition from 
the dipole configuration to the Saturn ring configuration 
occurring via this intermediate state. 




FIG. 12: Director map for the off-center ring configuration. 
The background contour plot represents the value of the order 
parameter. 
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FIG. 13: Dynamics of the transition from the satellite defect 
to the Saturn ring defect via an off-center ring configuration. 
We plot the z coordinate of the ring defect relative to the 
droplet center, versus number of the timesteps. As a guide, 
the dashed line corresponds to exponential decay with a time 
constant of 2 x 10 timesteps. 



IV. CONCLUSIONS 

We used molecular dynamics to study a spherical ob- 
ject, representing a small colloidal particle or droplet, 
suspended in a nematic liquid crystal matrix modelled 
with the soft ellipsoid potential. Homeotropic bound- 
ary conditions and strong anchoring create a hedgehog 
director configuration on the particle surface and in its 
vicinity; we have studied the defect structure around the 
particle which cancels this hedgehog defect and produces 
a uniform director field at large distances. We inves- 
tigated systems up to 10 6 particles with a domain de- 
composition program on a massively parallel computer, 
conducting runs of up to 1.5 x 10 6 timesteps. 

Our simulation results show that, in the case of a small 
droplet, only a Saturn ring defect is stable; an initial 



configuration containing a satellite defect evolves spon- 
taneously into a ring. For large enough droplets, how- 
ever, both the satellite and the ring defect are at least 
metastable on the timescale of the simulations. For inter- 
mediate sizes of colloid particle, a non-equatorial surface 
ring evolves very slowly towards the Saturn ring config- 
uration. 

Using order parameter maps we are able to resolve the 
position of the core of the defect; these also show that the 
nematic phase is strongly biaxial near the defect core. In 
these systems, the defect region is quite large: in some 
cases of a similar size to the colloid particle itself. Far 
from the particle the director profiles are in good agree- 
ment with those predicted by elastic theory. The elastic 
theory also predicts quite well the variation of defect po- 
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sition with colloid particle radius. 

It should be borne in mind that the colloid dimensions 
studied here are much smaller than those typically inves- 
tigated experimentally; indeed, it is quite surprising that 
the satellite defect is observable in our simulations, and 
that the off-center surface ring, which is never seen in 
the elastic theory, appears, albeit as a transition state. 
It is hoped that further work, aimed at estimating free 
energy differences between these structures, will clarify 
the situation. 
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